﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace NCM_MSTest.Alg
{
    public delegate double F(double x);
    public delegate double FF(double x1, double x2);

    /// <summary>
    /// 非线性方程（组）的数值解
    /// 对于超越方程和高次多项式的代数方程，通常使用逐次逼近的方法求根
    /// 二分法
    /// 迭代法
    /// </summary>
    internal class NumbericalSolutionOfNonlinearEquations
    {
        /// <summary>
        /// 确定有根区间
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="intervals"></param>
        public static void FindRootIntervals(F f, double a, double b
            , out double[][] intervals)
        {
            int n = 10;
            var temp = new List<Tuple<double, double>>();
            double h = (b - a) / n;
            double v = f(a);
            double xk = a + h;
            while (xk <= b)
            {
                double vk = f(xk);
                if (vk * v <= 0)
                    temp.Add(new Tuple<double, double>(xk - h, xk));
                xk += h;
                v = vk;
            }
            intervals = new double[temp.Count][];
            for (int i = 0; i < temp.Count; i++)
            {
                intervals[i] = new double[2] { temp[i].Item1, temp[i].Item2 };
            }
        }

        /// <summary>
        /// 二分法求根
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="root"></param>
        /// <param name="eps"></param>
        public static void Dichotomy(F f
            , double a, double b
            , out double root
            , double eps = 1e-5
            )
        {
            int n = 0;
            double r = f(a) * f(b);
            if (r == 0)
            {
                if (f(a) == 0)
                    root = a;
                else
                    root = b;
            }
            else if (r < 0)
            {
                n++;
                var t = (a + b) / 2;
                while (Math.Abs(f(t)) >= eps)
                {
                    if (f(a) * f(t) > 0)
                    {
                        a = t;
                    }
                    else
                    {
                        b = t;
                    }
                    n++;
                    t = (a + b) / 2;
                    Console.WriteLine(string.Format("n={0} f(a)={1} f(b)={2} t={3} f(t)={4}"
                        , n, f(a), f(b), t, f(t)));
                }
                root = t;
            }
            else
            {
                root = double.NaN;
                //无解
            }
        }

        /// <summary>
        /// 二分法求根 多根
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="roots"></param>
        /// <param name="eps"></param>
        public static void DichotomyEx(F f
           , double a, double b
           , out double[] roots
           , double eps = 1e-5
           )
        {
            double[][] intervals;
            FindRootIntervals(f, a, b, out intervals);
            if (intervals == null || intervals.Length == 0)
            {
                roots = null;
                return;
            }
            roots = new double[intervals.Length];
            for (int i = 0; i < intervals.Length; i++)
            {
                double da = intervals[i][0];
                double db = intervals[i][1];
                double r;
                Dichotomy(f, da, db, out r, eps);
                roots[i] = r;
            }
        }

        /// <summary>
        /// 不动点迭代法
        /// 将 y=f(x)改写成 x=g(x)
        /// 要求x*满足f(x*)=0,等价于求解x*使得x*=g(x*),则称x*为g(x)的不动点
        /// x_{n+1}=g(x_n),g(x)为迭代函数
        /// 不动点迭代收敛的必要条件是：|g'(x_0)| < 1 因此合理的迭代函数很重要
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="root"></param>
        /// <param name="eps"></param>
        public static void Iteration_FixedPoint(F f
            , double a, double b
            , int maxN //迭代的最大次数
            , out double root
            , double eps = 1e-5)
        {
            root = double.NaN;
            double t;
            F g;
            //确定迭代函数g(t)
            BuildIterationFunction(f, a, b, out g, out t);
            double t1 = g(t);
            int n = 0; //迭代次数
            //迭代
            while (Math.Abs(t1 - t) >= eps && n <= maxN)
            {
                t = t1;
                t1 = g(t);
                n++;
            }
            //输出t1 g(t1) n
            root = t1;
        }

        /// <summary>
        /// 不动点迭代法
        /// </summary>
        /// <param name="f"></param>
        /// <param name="g"></param>
        /// <param name="t"></param>
        /// <param name="maxN"></param>
        /// <param name="root"></param>
        /// <param name="n"></param>
        /// <param name="eps"></param>
        public static void Iteration_FixedPoint(F f
            , F g
            , double t
            , int maxN //迭代的最大次数
            , out double root
            , out int n
            , double eps = 1e-5)
        {
            root = double.NaN;
            //double t;
            //F g;
            ////确定迭代函数g(t)
            //BuildIterationFunction(f, a, b, out g, out t);
            double t1 = g(t);
            //int n = 0; //迭代次数
            n = 0;
            //迭代
            while (Math.Abs(t1 - t) >= eps && n <= maxN)
            {
                t = t1;
                t1 = g(t);
                n++;
            }
            //输出t1 g(t1) n
            root = t1;
        }

        /// <summary>
        /// 不动点迭代法
        /// </summary>
        /// <param name="g1"></param>
        /// <param name="g2"></param>
        /// <param name="tx"></param>
        /// <param name="maxN"></param>
        /// <param name="root"></param>
        /// <param name="n"></param>
        /// <param name="eps"></param>
        public static void Iteration_FixedPoint(
            FF g1
            , FF g2
            , double[] tx
            , int maxN //迭代的最大次数
            , out double[] root
            , out int n
            , double eps = 1e-4)
        {
            root = new double[2] { 0, 0 };
            double x01 = tx[0];
            double x02 = tx[1];
            //double t;
            //F g;
            ////确定迭代函数g(t)
            //BuildIterationFunction(f, a, b, out g, out t);
            double tx1 = g1(x01, x02);
            double tx2 = g2(x01, x02);
            //int n = 0; //迭代次数
            n = 1;
            Console.WriteLine($"{n}  {x01:F5}  {g1(x01, x02):F5} {x02:F5}  {g1(x01, x02):F5}");
            //迭代
            while ((Math.Abs(tx1 - x01) >= eps
                || Math.Abs(tx2 - x02) >= eps)
                && n <= maxN)
            {
                x01 = tx1;
                x02 = tx2;
                tx1 = g1(x01, x02);
                tx2 = g2(x01, x02);
                n++;
                Console.WriteLine($"{n}  {x01:F5}  {tx1:F5} {x02:F5}  {tx2:F5}");
            }
            //输出t1 g(t1) n
            root[0] = tx1;
            root[1] = tx2;
        }

        /// <summary>
        /// 构建迭代函数g(x)
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <returns></returns>
        public static bool BuildIterationFunction(F f, double a, double b, out F g, out double t)
        {
            t = (b + a) / 2;
            g = f;
            return true;
        }

        /// <summary>
        /// 牛顿迭代法
        /// 非线性方程f(x)在x0处有近似根，则在x0处的泰勒级数展开式：
        /// f(x)=f(x0)+f'(x0)(x-x0)+f''(x0)(x-x0)^2...
        /// 忽略高阶项：f(x)=f(x0)+f'(x0)(x-x0)
        /// 在点x0f附近：f(x0)+f'(x0)(x-x0)=0
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="maxN"></param>
        /// <param name="root"></param>
        /// <param name="eps"></param>
        public static void Iteration_Newton(F f
            , double t //自变量t的初始值
            , int maxN //迭代的最大次数
            , out double root
            , double eps = 1e-5)
        {
            root = double.NaN;
            //确定函数f、df
            F df;
            BuildDerivativeFunction(f, 1, out df);
            double t1 = t - f(t) / df(t);
            int n = 1;
            while (Math.Abs(f(t1)) >= eps && n <= maxN)
            {
                t = t1;
                t1 = t - f(t) / df(t);
                n++;
            }
            root = t1;
        }

        public static void Iteration_Newton(F f
            , F df
            , double t //自变量t的初始值
            , int maxN //迭代的最大次数
            , out double root
            , out int n
            , double eps = 1e-5)
        {
            root = double.NaN;
            //确定函数f、df
            //F df;
            //BuildDerivativeFunction(f, 1, out df);
            double t1 = t - f(t) / df(t);
            //int n = 1;
            n = 0;
            while (Math.Abs(f(t1)) >= eps && n <= maxN)
            {
                t = t1;
                t1 = t - f(t) / df(t);
                n++;
            }
            root = t1;
        }

        public static void Iteration_Newton(FF f1, FF f2
            , FF df11, FF df12
            , FF df21, FF df22
            , double[] t //自变量t的初始值
            , int maxN //迭代的最大次数
            , out double[] root
            , out int n
            , double eps = 1e-5)
        {
            root = new double[] { 0, 0 };
            double x01 = t[0];
            double x02 = t[1];
            n = 0;
            Console.WriteLine($"{n}  {x01:F5}  {f1(x01, x02):F5} {x02:F5}  {f2(x01, x02):F5}");
            double x11 = 0;
            double x12 = 0;
            while ((Math.Abs(f1(x01, x02)) >= eps
                || Math.Abs(f2(x01, x02)) >= eps)
                && n <= maxN)
            {
                n++;
                x11 = x01 -
                    (f1(x01, x02) * df22(x01, x02) - f2(x01, x02) * df12(x01, x02))
                    /
                    (df11(x01, x02) * df22(x01, x02) - df12(x01, x02) * df21(x01, x02));
                x12 = x02 -
                    (f2(x01, x02) * df11(x01, x02) - f1(x01, x02) * df21(x01, x02))
                    /
                    (df11(x01, x02) * df22(x01, x02) - df12(x01, x02) * df21(x01, x02));
                x01 = x11;
                x02 = x12;
                Console.WriteLine($"{n}  {x01:F5}  {x11:F5} {x02:F5}  {x12:F5}");
            }
            root[0] = x01;
            root[1] = x02;
        }

        /// <summary>
        /// 构建指定函数的指定阶数的导函数
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="order"></param>
        public static void BuildDerivativeFunction(F f, uint order
            , out F df)
        {
            df = f;
        }
    }
}
